## Daniel J. Mallinson
## Reproduction Code for "Who are Your Neighbors? The Role of Ideology and Decline of Geographic Proximity in the Diffusion of Policy Innovations" Published in Policy Studies Journal. 
## 2018-05-09
## 2019-01-09
## 2019-02-12

rm(list = ls(all = TRUE))

########################################################################
################## Load packages #######################################
########################################################################

install.packages(c("psych", "lme4", "reshape", "ggplot2", "gridExtra", "RColorBrewer", "scales", "MCMCpack"))

library(foreign)
library(psych)
library(lme4)
library(reshape2)
library(plyr)
library(ggplot2)
library(ggthemes)
library(scales)
library(splines)
library(gridExtra)
library(RColorBrewer)
library(MCMCpack)
library(Hmisc)
library(interplot)
library(lmtest)
library(sjPlot)

#Create color psublue
psublue <- rgb(red=0.0, green=0.14, blue=0.4, alpha=1)

########################################################################
################## Load data     #######################################
########################################################################

all.data <- read.csv("spid_macro_dataset_2019-01-10.csv")
all.data$X <- all.data$X.1 <- NULL
data <- all.data

#Add indicator for law/criminal justice policies
data$law <- 0
data$law[data$majortopic==12] <- 1

data$year_count <- data$year - 1960
data$time_log <- log(data$time+1)

#Keep data between 1960 and 2014
data <- data[which(data$year<=2014),]

##############################################################################
############# Select Variables and Create Decade Datasets ####################
##############################################################################

data[c("neighbor_prop","ideology_relative_hm","congress_majortopic", "percap_log","population_log",
                       "legprof_squire","mip","nyt", "time", "time2", "time3", "year_count")] <- scale(data[c("neighbor_prop","ideology_relative_hm","congress_majortopic", "percap_log","population_log",
                       "legprof_squire","mip","nyt", "time", "time2", "time3", "year_count")])


sixties <- data[which(data$first_year >= 1960 & data$first_year< 1970),]
seventies <- data[which(data$first_year>= 1970 & data$first_year< 1980),]
eighties <- data[which(data$first_year>=1980 & data$first_year< 1990), ]
nineties <- data[which(data$first_year>= 1990 & data$first_year< 2000),]
ots <- data[which(data$first_year>= 2000),]

##############################################################################
############################# Descriptives  ##################################
##############################################################################

describe(data, skew=FALSE, ranges=FALSE)
describe(sixties, skew=FALSE, ranges=FALSE)
describe(seventies, skew=FALSE, ranges=FALSE)
describe(eighties, skew=FALSE, ranges=FALSE)
describe(nineties, skew=FALSE, ranges=FALSE)
describe(ots, skew=FALSE, ranges=FALSE)

## Number of policies in overall analysis
length(unique(data$policy)) #556 policies

##############################################################################
########### Comparison to PA Agendas Project (SI) ############################
##############################################################################

## Pennsylvania Adoptions
data.pa <- read.csv("PA_Agendas_Act_Totals.csv")
head(data.pa)
data.pa[is.na(data.pa)] <- 0
head(data.pa)

total_col <- apply(data.pa[,-1], 1, sum)

pa.pcts <- lapply(data.pa[,-1], function(x){(x/total_col)*100})

pa.pcts<- as.data.frame(pa.pcts)
pa.pcts$year <- data.pa$Year

pa.pcts<- melt(pa.pcts, id.vars="year")
pa.pcts <- pa.pcts[which(pa.pcts$year<2015),]

labels <- c("Fiscal and Economic", "Civil Rights and Liberties", "Health", "Agriculture", "Labor and Employment",
		"Education", "Environment", "Energy", "Immigration", "Transportation", "Law", "Social Welfare",
		"Housing", "Finance", "Defense", "Science and Technology", "Foreign Trade", "International Affairs",
		"State Government Operations", "Public Lands", "Local Government")

png("pa_cumulative_plot.png", height=6, width=10, units="in", res=600)
ggplot(pa.pcts, aes(x=year, y=value, fill=variable))+
  geom_area(colour="black", size=.2, alpha=1)+
  scale_fill_manual(values=colorRampPalette(excel_new_pal()(6))(21), labels=labels) +
  labs(fill="Major Topics", x="Year", y="Percentage of All Topics")
dev.off()

## SPID Adoptions

adopts <- data[which(data$adopt==1),]
adopts <- subset(adopts, year >=1979 & year <=2014)
adopts <- adopts[c("majortopic", "year", "adopt")]

adopts.total <- as.data.frame.matrix(xtabs(adopt~year+majortopic,adopts))
year <- seq(1979,2014,1)
adopts.total <- cbind(year, adopts.total)
adopts.total$"24" <- 0
adopts.total.col <- apply(adopts.total[,-1], 1, sum)

adopts.pcts <- lapply(adopts.total[,-1], function(x){(x/adopts.total.col)*100})
adopts.pcts <- as.data.frame(adopts.pcts)

adopts.pcts$year <- year

adopts.pcts <- melt(adopts.pcts, id.vars="year")

png("spid_cumulative_plot.png", height=6, width=10, units="in", res=600)
ggplot(adopts.pcts, aes(x=year, y=value, fill=variable))+
  geom_area(colour="black", size=.2, alpha=1)+
  scale_fill_manual(values=colorRampPalette(excel_new_pal()(6))(21), labels=labels) +
  labs(fill="Major Topics", x="Year", y="Percentage of All Topics")
dev.off()

## Percentage of policies by major topic
pol.topics <- unique(data[c("policy", "majortopic")])
prop.table(table(pol.topics$majortopic))*100

##############################################################################
################# Functional Form Assumptions (SI) ###########################
##############################################################################

## Checking functional form assumption
model.logit <- glm(adopt ~ 
                        + neighbor_prop + ideology_relative_hm + congress_majortopic
                      + init_avail + init_qual + divided_gov + legprof_squire 
                      + percap_log + population_log
                      + mip + complexity_topic + mip*complexity_topic
                      + nyt + law + time, family=binomial(link="logit"), data=data)
summary(model.logit)

model.probit <- glm(adopt ~ 
                       + neighbor_prop + ideology_relative_hm + congress_majortopic
                      + init_avail + init_qual + divided_gov + legprof_squire 
                      + percap_log + population_log
                      + mip + complexity_topic + mip*complexity_topic
                      + nyt + law + time, family=binomial(link="probit"), data=data)

summary(model.probit)

lrtest(model.logit, model.probit)

model.cloglog <- glm(adopt ~ 
                        + neighbor_prop + ideology_relative_hm + congress_majortopic
                      + init_avail + init_qual + divided_gov + legprof_squire 
                      + percap_log + population_log
                      + mip + complexity_topic + mip*complexity_topic
                      + nyt + law + time, family=binomial(link="cloglog"), data=data)

summary(model.cloglog)

lrtest(model.logit, model.cloglog)

#N
nobs(model.logit) #299599

##############################################################################
################# Duration Dependence Assumptions (SI) #######################
##############################################################################

model.time <- glm(adopt ~ 
                        + neighbor_prop + ideology_relative_hm + congress_majortopic
                      + init_avail + init_qual + divided_gov + legprof_squire 
                      + percap_log + population_log
                      + mip + complexity_topic + mip*complexity_topic
                      + nyt + law + time, family=binomial(link="logit"), data=data)
summary(model.time)

model.logtime <- glm(adopt ~ 
                        + neighbor_prop + ideology_relative_hm + congress_majortopic
                      + init_avail + init_qual + divided_gov + legprof_squire 
                      + percap_log + population_log
                      + mip + complexity_topic + mip*complexity_topic
                      + nyt + law + time_log, family=binomial(link="logit"), data=data)

summary(model.logtime)

lrtest(model.time, model.logtime)

model.polytime <- glm(adopt ~ 
                        + neighbor_prop + ideology_relative_hm + congress_majortopic
                      + init_avail + init_qual + divided_gov + legprof_squire 
                      + percap_log + population_log
                      + mip + complexity_topic + mip*complexity_topic
                      + nyt + law + poly(time,3), family=binomial(link="logit"), data=data)
summary(model.polytime)

lrtest(model.time, model.polytime)

model.splinetime <- glm(adopt ~ 
                        + neighbor_prop + ideology_relative_hm + congress_majortopic
                      + init_avail + init_qual + divided_gov + legprof_squire 
                      + percap_log + population_log
                      + mip + complexity_topic + mip*complexity_topic
                      + nyt + law + ns(time, df=3), family=binomial(link="logit"), data=data)
summary(model.splinetime)

lrtest(model.time, model.splinetime)

model.nolaw <- glm(adopt ~ 
                        + neighbor_prop + ideology_relative_hm + congress_majortopic
                      + init_avail + init_qual + divided_gov + legprof_squire 
                      + percap_log + population_log
                      + mip + complexity_topic + mip*complexity_topic
                      + nyt + ns(time, df=3), family=binomial(link="logit"), data=data)
summary(model.nolaw)

lrtest(model.splinetime, model.nolaw)

##############################################################################
############################ Random Effects ##################################
##############################################################################

## Examine which random effect structure best fits the data using guidance from 
## Field, Andy, Jeremy Miles, and Zoe Field. 2012.  Discovering Statistics Using R. SAGE. 

interceptonly <- glm(adopt~1, family=binomial(link="logit"), data=data)
summary(interceptonly)

policyintercept <- glmer(adopt~ 1|policy, family=binomial(link="logit"), data=data, method="ML")
summary(policyintercept)

lrtest(interceptonly, policyintercept)

bothintercept <- glmer(adopt~ (1|policy) + (1|state), family=binomial(link="logit"), data=data, method="ML")
summary(bothintercept)

reint <- REsim(bothintercept)
plotREsim(reint)

lrtest(policyintercept, bothintercept)

randomneigh <- glmer(adopt~ neighbor_prop + (1+ neighbor_prop|policy), family=binomial(link="logit"), data=data)
summary(randomneigh)

lrtest(randomintercept, randomneigh)

randomboth <- glmer(adopt~ neighbor_prop + ideology_relative_hm + (1+ neighbor_prop + ideology_relative_hm |policy), family=binomial(link="logit"), data=data)
summary(randomboth)

lrtest(randomneigh, randomboth)

re <- REsim(randomboth)
plotREsim(re)

anova(policyintercept, interceptonly, randomneigh, randomboth) #random intercepts for policy are the better model

## Seems that adding random slopes for ideological distance improves model,
## but there is little variance in the random slopes

##############################################################################
############################ Table 1 #########################################
##############################################################################

## All Pooled Model
model.pooled <- glmer(adopt ~ 
                       neighbor_prop + ideology_relative_hm + congress_majortopic
                      + init_avail + init_qual + divided_gov + legprof_squire 
                      + percap_log + population_log
                      + mip + complexity_topic + mip*complexity_topic
                      + nyt + year_count + time_log + (1 + neighbor_prop | policy), 
                      control=glmerControl(optCtrl=list("optim", maxfun=5e5)), family=binomial(link="logit"), data=data)
summary(model.pooled)
ss <- getME(model.pooled,c("theta","fixef"))
model.pooled<- update(model.pooled,start=ss,control=glmerControl(optCtrl=list(maxfun=2e4)))
summary(model.pooled) 

se.pooled <- sqrt(diag(vcov(model.pooled)))
tab.pooled <- cbind(Est=fixef(model.pooled), LL=fixef(model.pooled) - 1.96*se.pooled, UL = fixef(model.pooled) + 1.96*se.pooled)
exp(tab.pooled)

## Year Interaction Model (Not reported in Table 1)
model.yearx <- glmer(adopt ~ 
                        neighbor_prop + ideology_relative_hm + 
				year_count + time_log + neighbor_prop*year_count + ideology_relative_hm*year_count +
				(1 + neighbor_prop| policy), 
                      control=glmerControl(optCtrl=list(maxfun=5e5)), family=binomial(link="logit"), data=data)
summary(model.yearx)
ss <- getME(model.yearx,c("theta","fixef"))
model.yearx<- update(model.yearx,start=ss,control=glmerControl(optCtrl=list(maxfun=2e4)))
summary(model.yearx)

se.yearx <- sqrt(diag(vcov(model.yearx)))
tab.yearx<- cbind(Est=fixef(model.yearx), LL=fixef(model.yearx) - 1.96*se.yearx, UL = fixef(model.yearx) + 1.96*se.yearx)
exp(tab.yearx)

## Interaction Full Model
model.full <- glmer(adopt ~ 
                   	neighbor_prop + ideology_relative_hm + congress_majortopic
                      + init_avail + init_qual + divided_gov + legprof_squire 
                      + percap_log + population_log
                      + mip + complexity_topic + mip*complexity_topic
                      + nyt + year_count + time_log + neighbor_prop*year_count + ideology_relative_hm*year_count
				              +(1 + neighbor_prop | policy), 
                      control=glmerControl(optCtrl=list(maxfun=5e5)), family=binomial(link="logit"), data=data)
summary(model.full)
ss <- getME(model.full,c("theta","fixef"))
model.full<- update(model.full,start=ss,control=glmerControl(optCtrl=list(maxfun=2e4)))
summary(model.full) #May need to repeat this several times for model to converge

se.full <- sqrt(diag(vcov(model.full)))
tab.full <- cbind(Est=fixef(model.full), LL=fixef(model.full) - 1.96*se.full, UL = fixef(model.full) + 1.96*se.full)
exp(tab.full)

##Plot of random effects
reEx <- REsim(model.full)
reEx$term[reEx$term=="ideology_relative_hm"] <- "Ideological Distance"
reEx$term[reEx$term=="neighbor_prop"] <- "Neighbor Adoptions"
reEx$groupFctr <- "Policy"

jpeg("fullmodel_re.jpg", width=8, height=6, units="in", res=600)
plotREsim(reEx)
dev.off()

##############################################################################
############################ Table 2 #########################################
##############################################################################

## Salience interaction model
model.salience <- glmer(adopt ~ 
                   	neighbor_prop + ideology_relative_hm + congress_majortopic
                      + init_avail + init_qual + divided_gov + legprof_squire 
                      + percap_log + population_log
                      + mip + complexity_topic + complexity_topic*mip
                      + nyt + year_count + time_log + ideology_relative_hm*mip
				+(1 + neighbor_prop | policy), 
                      control=glmerControl(optCtrl=list(maxfun=5e5)), family=binomial(link="logit"), data=data)
summary(model.salience)
ss <- getME(model.salience,c("theta","fixef"))
model.salience<- update(model.salience,start=ss,control=glmerControl(optCtrl=list(maxfun=2e4)))
summary(model.salience)

se.salience <- sqrt(diag(vcov(model.salience)))
tab.salience <- cbind(Est=fixef(model.salience), LL=fixef(model.salience) - 1.96*se.salience, UL = fixef(model.salience) + 1.96*se.salience)
exp(tab.salience)

## Complexity interaction model
model.complexity <- glmer(adopt ~ 
                   	neighbor_prop + ideology_relative_hm + congress_majortopic
                      + init_avail + init_qual + divided_gov + legprof_squire 
                      + percap_log + population_log
                      + mip + complexity_topic + complexity_topic*mip
                      + nyt + year_count + time_log + ideology_relative_hm*complexity_topic
				+(1 + neighbor_prop | policy), 
                      control=glmerControl(optCtrl=list(maxfun=5e5)), family=binomial(link="logit"), data=data)
summary(model.complexity)
ss <- getME(model.complexity,c("theta","fixef"))
model.complexity <- update(model.complexity,start=ss,control=glmerControl(optCtrl=list(maxfun=2e4)))
summary(model.complexity)

se.complexity <- sqrt(diag(vcov(model.complexity)))
tab.complexity <- cbind(Est=fixef(model.complexity), LL=fixef(model.complexity) - 1.96*se.complexity, UL = fixef(model.complexity) + 1.96*se.complexity)
exp(tab.complexity)

## Polarization Model
model.polarization <- glmer(adopt ~ 
                   	neighbor_prop + ideology_relative_hm + congress_majortopic
			    + avg_diffs
                      + init_avail + init_qual + divided_gov + legprof_squire 
                      + percap_log + population_log
                      + mip + complexity_topic + complexity_topic*mip
                      + nyt + year_count + time_log + ideology_relative_hm*avg_diffs + neighbor_prop*avg_diffs
				+(1 + neighbor_prop | policy), 
                      control=glmerControl(optCtrl=list(maxfun=5e5)), family=binomial(link="logit"), data=data)
summary(model.polarization)
ss <- getME(model.polarization,c("theta","fixef"))
model.polarization<- update(model.polarization,start=ss,control=glmerControl(optCtrl=list(maxfun=2e4)))
summary(model.polarization)

se.polarization <- sqrt(diag(vcov(model.polarization)))
tab.polarization <- cbind(Est=fixef(model.polarization), LL=fixef(model.polarization) - 1.96*se.polarization, UL = fixef(model.polarization) + 1.96*se.polarization)
exp(tab.polarization)

##############################################################################
################################ Figure 1 ####################################
##############################################################################

#https://cran.r-project.org/web/packages/interplot/vignettes/interplot-vignette.html 

#jpeg("figure1.jpg", height=8, width=11, units="in", res=600)
pdf("figure1.pdf", height=8, width=11)
a <- interplot(m=model.full, var1="complexity_topic", var2="mip", hist=TRUE, adjCI=TRUE) +
  xlab("Salience (MIP, Mean Centered)") + 
  ylab("Marginal Effect of Policy Complexity") + 
  theme_bw() + 
  geom_hline(yintercept=0, linetype="solid")

b <- plot_model(model.full, type="pred", terms=c("complexity_topic[0,1]", "mip[-0.669, 7.79]"), legend.title="Salience (MIP)",axis.title="Predicted Probabilities of Policy Adoption")
b[9]$labels$x <- "Complex Policy"
b[9]$labels$title <- ""
#model[9]$labels$y <- "Predicted Probability of Policy Adoption"

c <- interplot(m=model.full, var1="mip", var2="complexity_topic", hist=FALSE, adjCI=TRUE) +
  ylab("Marginal Effect of Salience (MIP)") + 
  #scale_x_discrete(name="Policy Complexity", labels=c("-0.7495224" = "Not Complex", "1.334102" = "Complex")) +
  scale_x_continuous(name="Policy Complexity", breaks=c(0, 1),
                     labels=c("Non-Complex", "Complex")) + 
  #xlim(-1, 1.5) +
  theme_bw() + 
  geom_hline(yintercept=0, linetype="solid")

d <- plot_model(model.full, type="pred", terms=c("mip", "complexity_topic[0,1]"), legend.title="Complex Policy",axis.title="Predicted Probabilities of Policy Adoption")
d[9]$labels$x <- "Salience (MIP, Mean Centered)"
d[9]$labels$title <- ""
#model[9]$labels$y <- "Predicted Probability of Policy Adoption"

grid.arrange(a,b,c,d,nrow=2)
dev.off()


##############################################################################
############################ Figure  2 #######################################
##############################################################################

#jpeg("figure2.jpg", height=8, width=11, units="in", res=600)
pdf("figure2.pdf", height=8, width=11)
e <- interplot(m=model.full, var1="neighbor_prop", var2="year_count", hist=TRUE, adjCI=TRUE) +
  xlab("Years, 1960 to 2014 (Mean Centered)") + 
  ylab("Marginal Effect of Neighbor Adoptions") + 
  theme_bw() + 
  geom_hline(yintercept=0, linetype="solid")

f <- plot_model(model.full, type="pred", terms=c("neighbor_prop", "year_count[-2.59, 1.93]"), legend.title="Year (Centered)",axis.title="Predicted Probabilities of Policy Adoption")
f[9]$labels$x <- "Proportion of Neighbors Adopting (Mean Centered)"
f[9]$labels$title <- ""
#model[9]$labels$y <- "Predicted Probability of Policy Adoption"

g <- interplot(m=model.full, var1="ideology_relative_hm", var2="year_count", hist=FALSE, adjCI=TRUE) +
  ylab("Marginal Effect of Ideological Distance") + 
  xlab("Years, 1960 to 2014 (Mean Centered)") + 
  theme_bw() + 
  geom_hline(yintercept=0, linetype="solid")

h <- plot_model(model.full, type="pred", terms=c("ideology_relative_hm", "year_count[-2.59, 1.93]"), legend.title="Year (Centered)",axis.title="Predicted Probabilities of Policy Adoption")
h[9]$labels$x <- "Ideological Distance (Mean Centered)"
h[9]$labels$title <- ""
#model[9]$labels$y <- "Predicted Probability of Policy Adoption"

grid.arrange(e,f,g,h,nrow=2)
dev.off()

##############################################################################
############################ Figure  3 #######################################
##############################################################################

#jpeg("figure3.jpg", height=8, width=11, units="in", res=600)
pdf("figure3.pdf", height=8, width=11)
i <- interplot(m=model.salience, var1="ideology_relative_hm", var2="mip", hist=TRUE, adjCI=TRUE) +
  xlab("Salience (MIP, Mean Centered)") + 
  ylab("Marginal Effect of Ideological Distance") + 
  theme_bw() + 
  geom_hline(yintercept=0, linetype="solid")

j <- plot_model(model.salience, type="pred", terms=c("ideology_relative_hm", "mip[-0.669, 7.79]"), legend.title="MIP (Centered)",axis.title="Predicted Probabilities of Policy Adoption")
j[9]$labels$x <- "Ideological Distance (Mean Centered)"
j[9]$labels$title <- ""
#model[9]$labels$y <- "Predicted Probability of Policy Adoption"

k <- interplot(m=model.complexity, var1="ideology_relative_hm", var2="complexity_topic", hist=FALSE, adjCI=TRUE) +
  ylab("Marginal Effect of Ideological Distance") + 
  scale_x_continuous(name="Policy Complexity", breaks=c(0, 1),
                     labels=c("Non-Complex", "Complex")) + 
  theme_bw() + 
  geom_hline(yintercept=0, linetype="solid")

l <- plot_model(model.full, type="pred", terms=c("ideology_relative_hm", "complexity_topic[0,1]"), legend.title="Complex Policy",axis.title="Predicted Probabilities of Policy Adoption")
l[9]$labels$x <- "Ideological Distance (Mean Centered)"
l[9]$labels$title <- ""
#model[9]$labels$y <- "Predicted Probability of Policy Adoption"

grid.arrange(i,j,k,l,nrow=2)
dev.off()

##############################################################################
############################ Figure  4 #######################################
##############################################################################

#jpeg("figure4.jpg", height=8, width=11, units="in", res=600)
pdf("figure4.pdf", height=8, width=11)
m <- interplot(m=model.polarization, var1="neighbor_prop", var2="avg_diffs", hist=TRUE, adjCI=TRUE) +
  xlab("Legislative Polarization (Mean Centered)") + 
  ylab("Marginal Effect of Neighbor Adoptions") + 
  theme_bw() + 
  geom_hline(yintercept=0, linetype="solid")

n <- plot_model(model.polarization, type="pred", terms=c("neighbor_prop", "avg_diffs[0.2, 3.165]"), legend.title="Polarization",axis.title="Predicted Probabilities of Policy Adoption")
n[9]$labels$x <- "Neighbor Adoptions (Mean Centered)"
n[9]$labels$title <- ""
#model[9]$labels$y <- "Predicted Probability of Policy Adoption"

o <- interplot(m=model.polarization, var1="ideology_relative_hm", var2="avg_diffs", hist=FALSE, adjCI=TRUE) +
  xlab("Legislative Polarization (Mean Centered)") + 
  ylab("Marginal Effect of Ideological Distance") + 
  theme_bw() + 
  geom_hline(yintercept=0, linetype="solid")

p <- plot_model(model.polarization, type="pred", terms=c("ideology_relative_hm", "avg_diffs[0.2, 3.165]"), legend.title="Polarization", axis.title="Predicted Probabilities of Policy Adoption")
p[9]$labels$x <- "Ideological Distance (Mean Centered)"
p[9]$labels$title <- ""
#model[9]$labels$y <- "Predicted Probability of Policy Adoption"

grid.arrange(m,n,o,p,nrow=2)
dev.off()

##############################################################################
############################ Figure 5 ########################################
##############################################################################

#https://www.rdocumentation.org/packages/lme4/versions/1.1-17/topics/ranef
#https://stackoverflow.com/questions/13847936/in-r-plotting-random-effects-from-lmer-lme4-package-using-qqmath-or-dotplot
#http://www.strengejacke.de/sjPlot/reference/sjp.lmer.html
#https://rdrr.io/github/davebraze/FDB1/man/ggCaterpillar.html
#https://stats.idre.ucla.edu/r/dae/mixed-effects-logistic-regression/
#https://stats.idre.ucla.edu/other/mult-pkg/introduction-to-generalized-linear-mixed-models/

randoms <- ranef(model.full, condVar=TRUE)
pv <- attr(randoms[[1]],"postVar")
vc <- split(pv, slice.index(pv,3))
vcdata <- as.data.frame(vc)
vcdata <- vcdata[c(1,4),]
vcdata <- t(vcdata)
vcdata <- sqrt(vcdata)
vcdata <- as.data.frame(vcdata)
colnames(vcdata) <- c("intercept", "neighbor_prop")

coeffs.neighbor <- data.frame(policy=rownames(coef(model.full)$policy),
			coefficients=coef(model.full)$policy["neighbor_prop"],
			error = 0.02269+vcdata$neighbor_prop,
			pos = as.vector(seq(1,1112,2)))
coeffs.neighbor$uci <- coeffs.neighbor$neighbor_prop + 1.96*coeffs.neighbor$error
coeffs.neighbor$lci <- coeffs.neighbor$neighbor_prop - 1.96*coeffs.neighbor$error
for(i in 1:556){
if(sign(coeffs.neighbor$uci[i])==sign(coeffs.neighbor$lci[i])){
	coeffs.neighbor$sig[i] = 1}
	else{
	coeffs.neighbor$sig[i] = 0}}
		
sig.neigh <- coeffs.neighbor[which(coeffs.neighbor$sig==1),]
sig.neigh <- sig.neigh[order(-sig.neigh$neighbor_prop),]
length <- nrow(sig.neigh)


jpeg("neighbor_effects1.jpg", height=8, width=5, units="in", res=600)
par(mar=c(5,9,3,1))
plot(x=sig.neigh$neighbor_prop[1:74], y=seq(1,74,1), xlim=c(-2,2), axes=FALSE, pch=20, ylab="", xlab="Policy-Level Effect", main="Neighbor Adoptions")
segments(sig.neigh$lci[1:74], seq(1,74,1), sig.neigh$uci[1:74], seq(1,74,1))
abline(v=0)
axis(1, at=seq(-2,2,.5), labels=seq(-2,2,.5))
labels <- c("Amber Alert", "Zero Tolerance", "Renewable Energy Standards", "Grandparent Visitation", "DUI", 
		"Underage DUI .02", "Smoking Ban (Bars)", "Any Photo ID", "Meth Precursor Ban", "Smoking Ban (Restaurant)",
		"Film Tax Credit", "Medical Abortion Restriction", "Doctor Gag Ban", "Megans Law", "College Savings",
		"College Savings (2)", "Enviro Building Standards", "Medical Savings Accounts", "UCC Article 9", "Autism Coverage",
		"Civil Unions", "Paper Terrorism", "Drug Ban (Restaurant)", "Anti-Stalking", "Statewide Smoking Ban", 
		"Zero Tolerance", "Research on Ru486", "Minimum Wage Above Federal", "Motor Voter", "Human Trafficking", 
		"UCC Article 9 Amendments", "Public Benefits", "Discrimination (Sexual Orientation)", "Electric Deregulation",
		"Fetal Ultrasound", "Gender Discrimination EO", "Performance-Based Funding", "Emergency Contraception", 
		"Anti-Bullying", "Car Emissions", "Abortion Public Funding Ban", "Credit Freeze", "Medical Marijuana", 
		"Good Samaritan", "Multistate Highway Agreement", "Absentee Voting", "Indian Gaming", "Breast Density Notification", 
		"Individual Development Accounts", "Lemon Law", "Dependent Coverage Expansion", "Electronic Waste",
		"Hate Crime (LGBT)", "Renewable Portfolio Standards", "Any Voter ID", "DWI Reform", "Pain Management", 
		"Employment Discrimination (LGBT)", "Child Car Restraint", "Fetal Remains Disposal", "Keg Registration", 
		"Limited Liability Company", "Syringe Exchange", "Fetal Homicide", "Bike Helmets", "Public Benefit Fund", 
		"Corporal Punishment", "Gay Marriage", "Net Metering", "Electric Restructuring", "Performance-Based Funding",
		"Physician Abortion Refusal", "Instate Tuition (Immigrant)", "Interstate Discovery")
axis(2, at=seq(1,74,1), labels=labels, las=2, cex.axis=.6)
dev.off()

jpeg("neighbor_effects2.jpg", height=8, width=5, units="in", res=600)
par(mar=c(5,9,3,1))
plot(x=sig.neigh$neighbor_prop[75:148], y=seq(1,74,1), xlim=c(-2,2), axes=FALSE, pch=20, ylab="", xlab="Policy-Level Effect", main="Neighbor Adoptions")
segments(sig.neigh$lci[75:148], seq(1,74,1), sig.neigh$uci[75:148], seq(1,74,1))
abline(v=0)
axis(1, at=seq(-2,2,.5), labels=seq(-2,2,.5))
labels <- c("Victim Notification Req", "Anti-Hazing", "State Victim Notification System", "Dam Removal", "Concealed Carry",
		"Charter Schools", "Terrorism Funding", "Living Wills", "Strategic Planning (Transportation)", 
		"Stand Your Ground", "Victim Rights", "Identity Theft", "Lottery", "DUI08", "Employment Discrimination (LBGT)", 
		"Infant Hearing Screen", "Imitation Drug Regulation", "Public Accomidation", "Interstate Family Support (2001)",
		"Interstate Family Support (1992)", "Eminent Scholar", "English Language", "Post-Convication DNA",
		"Insurance Regulation Compact", "Bike Helmet", "Anti-Stalking", "DNA for Exoneration", "Smoking Ban (Workplaces)",
		"Prudent Management Act", "Nonresident Violator Compact", "Any Voter ID", "Doctor Apology", "Breastfeeding in Public",
		"Zero Tolerance", "Repeal Anti-Sodomy", "UCC Article 8", "Child Abuse", "Son of Sam", "Net Metering",
		"TOD Security Registration", "Smoking Ban (Workplace)", "Victim's Rights Amendment", "Prudent Investor Act",
		"UCC Article 1", "Collective Bargaining", "Colan Cancer Screening", "Preexisting Conditions",
		"Open Carry", "WMD Ban", "Public Breast Feeding",
		"Pro-Life License Plates", "Graduated Driver's License", "Licensed Physician for Abortion", "Victim's Compensation",
		"Lottery", "Son of Sam", "UCC Article 2a", "Renewal of Health Insurance", "Health Insurance Portability", 
		"R&D Tax Credit", "EMA Compact", "Hate Crimes", "Dual Enrollment", "University Record System", "Discrimination Ban (Disabilities)",
		"Principal and Income Act", "Hate Crime (Minorities)", "Placement of Children Compact", "Determination of Death", 
		"Community Service Mandate", "Post-Viability Abortion Ban", "Motorcyle Helmets", "Fair Employment (Race)",
		"Anatomical Gift")
axis(2, at=seq(1,74,1), labels=labels, las=2, cex.axis=.6)
dev.off()

jpeg("figure5.jpg", height=10, width=8, units="in", res=600)
par(mfrow=c(1,2))
par(mar=c(5,9,3,1))
plot(x=sig.neigh$neighbor_prop[75:148], y=seq(1,74,1), xlim=c(-2,2), axes=FALSE, pch=20, ylab="", xlab="Policy-Level Effect", main="Neighbor Adoptions")
segments(sig.neigh$lci[75:148], seq(1,74,1), sig.neigh$uci[75:148], seq(1,74,1))
abline(v=0)
axis(1, at=seq(-2,2,.5), labels=seq(-2,2,.5))
labels <- c("Victim Notification Req", "Anti-Hazing", "State Victim Notification System", "Dam Removal", "Concealed Carry",
		"Charter Schools", "Terrorism Funding", "Living Wills", "Strategic Planning (Transportation)", 
		"Stand Your Ground", "Victim Rights", "Identity Theft", "Lottery", "DUI08", "Employment Discrimination (LBGT)", 
		"Infant Hearing Screen", "Imitation Drug Regulation", "Public Accomidation", "Interstate Family Support (2001)",
		"Interstate Family Support (1992)", "Eminent Scholar", "English Language", "Post-Convication DNA",
		"Insurance Regulation Compact", "Bike Helmet", "Anti-Stalking", "DNA for Exoneration", "Smoking Ban (Workplaces)",
		"Prudent Management Act", "Nonresident Violator Compact", "Any Voter ID", "Doctor Apology", "Breastfeeding in Public",
		"Zero Tolerance", "Repeal Anti-Sodomy", "UCC Article 8", "Child Abuse", "Son of Sam", "Net Metering",
		"TOD Security Registration", "Smoking Ban (Workplace)", "Victim's Rights Amendment", "Prudent Investor Act",
		"UCC Article 1", "Collective Bargaining", "Colan Cancer Screening", "Preexisting Conditions",
		"Open Carry", "WMD Ban", "Public Breast Feeding",
		"Pro-Life License Plates", "Graduated Driver's License", "Licensed Physician for Abortion", "Victim's Compensation",
		"Lottery", "Son of Sam", "UCC Article 2a", "Renewal of Health Insurance", "Health Insurance Portability", 
		"R&D Tax Credit", "EMA Compact", "Hate Crimes", "Dual Enrollment", "University Record System", "Discrimination Ban (Disabilities)",
		"Principal and Income Act", "Hate Crime (Minorities)", "Placement of Children Compact", "Determination of Death", 
		"Community Service Mandate", "Post-Viability Abortion Ban", "Motorcyle Helmets", "Fair Employment (Race)",
		"Anatomical Gift")
axis(2, at=seq(1,74,1), labels=labels, las=2, cex.axis=.6)

par(mar=c(5,9,3,1))
plot(x=sig.neigh$neighbor_prop[1:74], y=seq(1,74,1), xlim=c(-2,2), axes=FALSE, pch=20, ylab="", xlab="Policy-Level Effect", main="Neighbor Adoptions")
segments(sig.neigh$lci[1:74], seq(1,74,1), sig.neigh$uci[1:74], seq(1,74,1))
abline(v=0)
axis(1, at=seq(-2,2,.5), labels=seq(-2,2,.5))
labels <- c("Amber Alert", "Zero Tolerance", "Renewable Energy Standards", "Grandparent Visitation", "DUI", 
		"Underage DUI .02", "Smoking Ban (Bars)", "Any Photo ID", "Meth Precursor Ban", "Smoking Ban (Restaurant)",
		"Film Tax Credit", "Medical Abortion Restriction", "Doctor Gag Ban", "Megans Law", "College Savings",
		"College Savings (2)", "Enviro Building Standards", "Medical Savings Accounts", "UCC Article 9", "Autism Coverage",
		"Civil Unions", "Paper Terrorism", "Drug Ban (Restaurant)", "Anti-Stalking", "Statewide Smoking Ban", 
		"Zero Tolerance", "Research on Ru486", "Minimum Wage Above Federal", "Motor Voter", "Human Trafficking", 
		"UCC Article 9 Amendments", "Public Benefits", "Discrimination (Sexual Orientation)", "Electric Deregulation",
		"Fetal Ultrasound", "Gender Discrimination EO", "Performance-Based Funding", "Emergency Contraception", 
		"Anti-Bullying", "Car Emissions", "Abortion Public Funding Ban", "Credit Freeze", "Medical Marijuana", 
		"Good Samaritan", "Multistate Highway Agreement", "Absentee Voting", "Indian Gaming", "Breast Density Notification", 
		"Individual Development Accounts", "Lemon Law", "Dependent Coverage Expansion", "Electronic Waste",
		"Hate Crime (LGBT)", "Renewable Portfolio Standards", "Any Voter ID", "DWI Reform", "Pain Management", 
		"Employment Discrimination (LGBT)", "Child Car Restraint", "Fetal Remains Disposal", "Keg Registration", 
		"Limited Liability Company", "Syringe Exchange", "Fetal Homicide", "Bike Helmets", "Public Benefit Fund", 
		"Corporal Punishment", "Gay Marriage", "Net Metering", "Electric Restructuring", "Performance-Based Funding",
		"Physician Abortion Refusal", "Instate Tuition (Immigrant)", "Interstate Discovery")
axis(2, at=seq(1,74,1), labels=labels, las=2, cex.axis=.6)
dev.off()

##############################################################################
######################### SI Figures 3a-d ####################################
##############################################################################

## Plot of non-significant effects for SI
insig.neigh <- coeffs.neighbor[which(coeffs.neighbor$sig==0),]
insig.neigh <- insig.neigh[order(-insig.neigh$neighbor_prop),]
length <- nrow(insig.neigh)


jpeg("si_fig3a.jpg", height=10, width=8, units="in", res=600)
par(mar=c(5,20,3,1))
plot(x=insig.neigh$neighbor_prop[1:102], y=seq(1,102,1), xlim=c(-2,2), axes=FALSE, pch=20, ylab="", xlab="Policy-Level Effect", main="Neighbor Adoptions")
segments(insig.neigh$lci[1:102], seq(1,102,1), insig.neigh$uci[1:102], seq(1,102,1))
abline(v=0)
axis(1, at=seq(-2,2,.5), labels=seq(-2,2,.5))
axis(2, at=seq(1,102,1), labels=insig.neigh$policy[1:102], las=2, cex.axis=.6)
dev.off()

jpeg("si_fig3b.jpg", height=10, width=8, units="in", res=600)
par(mar=c(5,17,3,1))
plot(x=insig.neigh$neighbor_prop[103:204], y=seq(1,102,1), xlim=c(-2,2), axes=FALSE, pch=20, ylab="", xlab="Policy-Level Effect", main="Neighbor Adoptions")
segments(insig.neigh$lci[103:204], seq(1,102,1), insig.neigh$uci[103:204], seq(1,102,1))
abline(v=0)
axis(1, at=seq(-2,2,.5), labels=seq(-2,2,.5))
axis(2, at=seq(1,102,1), labels=insig.neigh$policy[103:204], las=2, cex.axis=.6)
dev.off()

jpeg("si_fig3c.jpg", height=10, width=8, units="in", res=600)
par(mar=c(5,20,3,1))
plot(x=insig.neigh$neighbor_prop[205:306], y=seq(1,102,1), xlim=c(-2,2), axes=FALSE, pch=20, ylab="", xlab="Policy-Level Effect", main="Neighbor Adoptions")
segments(insig.neigh$lci[205:306], seq(1,102,1), insig.neigh$uci[205:306], seq(1,102,1))
abline(v=0)
axis(1, at=seq(-2,2,.5), labels=seq(-2,2,.5))
axis(2, at=seq(1,102,1), labels=insig.neigh$policy[205:306], las=2, cex.axis=.6)
dev.off()

jpeg("si_fig3d.jpg", height=10, width=8, units="in", res=600)
par(mar=c(5,20,3,1))
plot(x=insig.neigh$neighbor_prop[307:408], y=seq(1,102,1), xlim=c(-2,2), axes=FALSE, pch=20, ylab="", xlab="Policy-Level Effect", main="Neighbor Adoptions")
segments(insig.neigh$lci[307:408], seq(1,102,1), insig.neigh$uci[307:406], seq(1,102,1))
abline(v=0)
axis(1, at=seq(-2,2,.5), labels=seq(-2,2,.5))
axis(2, at=seq(1,102,1), labels=insig.neigh$policy[307:408], las=2, cex.axis=.6)
dev.off()


##############################################################################
############################# Figure (Unused)  ###############################
##############################################################################

## Stacked bar chart of policies by decade
# Pull adoptions for each decade
sixties.adopt <- sixties[which(sixties$adopt==1),]
seventies.adopt <- seventies[which(seventies$adopt==1),]
eighties.adopt <- eighties[which(eighties$adopt==1),]
nineties.adopt <- nineties[which(nineties$adopt==1),]
ots.adopt <- ots[which(ots$adopt==1),]

# Include policy name, state, year, adopt, and major topic code
barchart.vars <- c("adopt", "majortopic")
sixties.adopt <- sixties.adopt[barchart.vars]
sixties.adopt$decade <- 1960
nrow(sixties.adopt)

seventies.adopt <- seventies.adopt[barchart.vars]
seventies.adopt$decade <- 1970
nrow(seventies.adopt)

eighties.adopt <- eighties.adopt[barchart.vars]
eighties.adopt$decade <- 1980
nrow(eighties.adopt)

nineties.adopt <- nineties.adopt[barchart.vars]
nineties.adopt$decade <- 1990
nrow(nineties.adopt)

ots.adopt <- ots.adopt[barchart.vars]
ots.adopt$decade <- 2000
nrow(ots.adopt)

#Create single dataset of adoptions
barchart.data <- rbind(sixties.adopt, seventies.adopt, eighties.adopt, nineties.adopt, ots.adopt)

#Calculate proportions by decade

test <- aggregate(barchart.data, by=list(barchart.data$decade, barchart.data$majortopic), FUN=sum)
test <- test[c("Group.1", "Group.2", "adopt")]
test2 <- ddply(test, "Group.1", transform, prop=adopt/sum(adopt)*100)

# create a ggplot  

#pdf("time_topics.pdf", height=8.5, width=11)
jpeg("time_topics.jpg", height=8.5, width=11, units="in", res=1000)
ggplot(test2, aes(x=as.factor(test2$Group.1), y=test2$prop, fill=as.factor(test2$Group.2))) + 
  # plot the bars that are black
  #geom_bar(stat="identity", colour="black") + 
  geom_bar(position="fill", stat="identity") +
  #Change legend font size, title and labels
  scale_fill_manual(name="Major Topics",
				values=colorRampPalette(brewer.pal(12,"Paired"))(20),
				breaks=as.character(c(1,2,3,5,6,7,8,10,12,13,14,15,20,21,24)),
				labels=c("Fiscal", "Civil Rights", "Health", "Labor", "Education",
				"Environment", "Energy", "Transportation", "Law and Crime",
				"Social Welfare", "Housing", "Commerce", "Administrative", 
				"Public Resource Mgmt", "Local Government")) +
  # reverse the legend so it matches the colors
  guides(fill=guide_legend(reverse=TRUE)) + 
  #Gray fill behind legend
  theme(legend.background = element_rect(fill="gray90", size=.5, linetype="dotted"), legend.text=element_text(size=24)) +
  #Change x and y axis labels
  scale_x_discrete("Decade", breaks=c("1960", "1970", "1980", "1990", "2000"), labels=c("1960s \n (n=166)", "1970s \n(n=364)", "1980s \n (n=735)", "1990s \n (n=1343)", "2000s \n (n=780)")) +
  scale_y_continuous("Percentage of Total Adoptions within Decade", labels=percent)+
  theme_bw()+
  theme(
    plot.background = element_blank()
   ,panel.grid.major = element_blank()
   ,panel.grid.minor = element_blank()
   ,panel.border = element_blank()
  )+
  theme(axis.line = element_line(color = 'black'))+
  theme(axis.text=element_text(size=12),
        axis.title=element_text(size=14,face="bold"))
dev.off() 

# See for more info: http://www.cookbook-r.com/Graphs/Legends_%28ggplot2%29/#changing-the-position-of-the-legend
# Quick color reference: http://sape.inf.usi.ch/quick-reference/ggplot2/colour
# Generating colors using RColorBrewer http://www.r-bloggers.com/r-using-rcolorbrewer-to-colour-your-figures-in-r/

